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Problem Definition 


Program COORD has been developed to generate a ’C'-type, two-dimensional 
coordinate system for isolated airfoil calculations. It is to be used in 
conjunction with a Navier— Stokes or Guler analysis. The coordinate system is 
an important component of. such analyses. An inappropriate coordinate system 
may lead to difficulty in obtaining a converged solution and even exhibit 
physically unrealistic predictions due to geometric truncation error. 

Therefore, generation of a viable system is mandatory. Any coordinate 
system used in the analysis should satisfy conditions of (i) generality, 

(ii) smoothness, (iii) resolvability and (iv) allow easy application of 
boundary conditions. Obviously, a coordinate system must be sufficiently 
general to allow application to the problem of interest if it is to be 
practical. The metric data associated with the coordinate system must be 
sufficiently smooth so that the variation from grid point to grid point 
does not lead to a numerical solution dominated by metric coefficient 
truncation error; it should be noted that this requirement differs from the 
requirement of the existence of a specified number of transformation 
derivatives. The coordinate system must resolve flow regions where rapid flow 
field changes occur. Finally, coordinates should allow accurate implementation 
of boundary conditions. The coordinate system used here satisfies all these 
requirements. 

In developing coordinate system generation methods several possibilities 
exist. These are approaches based upon conformal methods, approaches based 
upon solutions of a Poisson equation such as those of Ref. 1 and constructive 
approaches such as that of Refs. 2 and 3. The present procedure is a 
constructive approach, and is described in the following section. 

Method of Solution 

The coordinate system consists of a set of two families of curves; the 
£ — constant curves such as lines G'G in Fig. 1, and the q = constant curves 
such as ABCD or A'ED' in Fig. 1. Under the present effort the coordinate 
generation process has been automated and the coordinate construction process 
is as follows. The coordinate system is constructed by first forming the inner 
loop A'ED' which includes the blade. The blade shape is specified either via 
an equation in SUBROUTINE PRFL or via discrete coordinate locations. Since the 
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code constructs a f C f -type grid, the actual trailing edge region is replaced by 
a cusped trailing edge. This modification in geometry is not expected to 
significantly affect most airfoil flow field calculations for two reasons. 
First, many airfoils have relatively sharp trailing edges, and in these cases 
the cusped approximation represents a small approximation to the actual 
geometry. Secondly, when the coordinate system is used with a Navier-Stokes 
analysis, boundary layers appear on the blade and often the boundary layer 
displacement effects are much larger than the airfoil thickness itself in the 
trailing edge region. In these cases the cusped geometry approximation has 
little effect on the flow field calculations. However, there are some 
instances in which the cusped trailing edge may make a significant difference. 
One such case occurs when the near wake region is the region of high interest. 
In such cases the f C f grid cusped trailing edge approximation should be used 
with caution. 

Construction of the inner loop is followed by constructing an outer loop 
AGBCD. Both the inner and outer loops are then represented by parametric 
curves 


x = x(s) , y = y(s) 


( 1 ) 


where the parameter varies from zero to unity. The present coordinate 
generation process utilizes a multi-part transformation for the inner loop. 
First x and y are expressed as a function of s', the physical distance along 
the curve. After normalizing s' such that its range is between zero and unity, 
a local parametrization focuses upon the leading edge region and places 
s" = 0.5 at the leading edge point where highest resolution is sought. 

This location must be specified by the user and in general will be in the 
vicinity of the expected front stagnation point. This initial transformation 
is done automatically once the location where maximum resolution is required is 
specified. Following this first parametrization, s is related to s'* via a 
hyperbolic tangent parameterization centered about the leading edge point where 
s" = 0.5. 

In this transformation for s" <0.5 


s 


I I + ks" 
2T " I - ks" 


(2a) 
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and for s" >0.5 


* * I - 



In 


I + k(l-s") 

7 — k ( I — a M ) 


(2b) 


where 


k.2tanh(I.) ( 2c) 

and Che paramecer T controls the grid distribution. Increasing T packs points 
in the vicinity of s" = 0.5. 

The outer loop AGBCD is then parameterized via a cubic transformation in 
which coefficients of the cubic are related to the inner loop transformation. 
Two additional transformations follow. First an exponential stretch is imposed 
on both the inner and outer loops downstream of the trailing edge. This 
transformation is of the form 

s ' = S L |n [! + G(S L - s)] S<s L ( 3 ) 

where S L is the parameter value at the beginning of the transformation on 
both the inner and outer loop. Similarly, the parameter is stretched in an 
analogous manner for s > (1 - S L ) so that for corresponding points M and N on 
opposite sides of the outer loop 


where T and L are on opposite sides of the branch cut or outer loop as shown in 
Fig. 1. Finally, a local transformation is performed to pack points in the 

region of the trailing edge on the inner loop and in the corresponding region 
in the outer loop. 

Following construction and parameterization of the inner and outer primary 
loops, two additional loops are formed which are termed the inner and outer 
intermediate or secondary loops. These loops are constructed as shown in 
Fig. 2. The first intermediate loop is constructed downstream of x/c = 0.04 a 
given vertical distance, hi, from the inner loop. Upstream of x/c = 0.04, the 
distance is measured along a pseudo-radial line with center at x/c = .04 on the 
camber line as shown in Fig. 2. The outer intermediate loop is constructed 
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in an analogous manner at a distance V12 from the outer loop. The inner and 
outer secondary loops are parameterized to correspond to the primary loop 
points they are connected with via the vertical or pseudo-radial lines; i.e. , 


Similarly, 


S B' * S B 


S T'= S T 


These four parameterized loops allow the construction of the pseudo-radial 
lines such as G’G of Fig. 1 via the multi-loop method originally developed by 
Eiseman (Ref. 4). 

The multi-loop method requires introduction of a position vector P(r,s) 
with components (x, y) which will represent the coordinate location of grid 
points. Based upon the four-loop construction process, vectors Pi(s) are 
defined with i = 1,2, 3, 4. Each Pf has coordinate (x^, yj) associated 
with it at specific values of s through the previously described construction 
process. A radial paramter, r, is then introduced. This paramter is defined 
as the distance from the loop in question to the inner loop normalized by the 
distance from the outer loop to the inner loop. Thus, ri = 0, r2 “ hi/h3, 
r 3 " 0*3 “ b2)/h3, r4 = 1. With the definition of these quantities the general 
position vector P(r,s) is related to the loop position vectors Pi(s), P2(s), 
P3<s) and ?4(s) via 


P(r,s) «(l-r) 2 (l-a,r)^(s) + (a,+ 2)(l-r) 2 rP 2 (s) 

+ r 2 [l-a 2 (l-r)]^(s) + (a 2 +2)r 2 (l-r)P 3 (s) 


(4) 


where 



2 

a 2 = 3(l-r 2 )-| 


(5) 
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It should be noted that at r = 0, P(0,s) = Pi(s) and at r = 1, P(l,s) = P 4 (s). 
Further since at r 8 0, 

*^jr (0,s) - [p 2 (s)-p ( (s)](a,+2) (6) 

and at r = 1 

■^•(l,s) 8 [^(s)-^(s)](a z +2) ( 7 ) 

specification of the derivatives at the inner and outer boundaries is 
determined by the parametric representation of intermediate loops 2 and 3. 

Thus, the four-loop method allows specification of the boundary point locations 
and coordinate angles at these boundaries. 

After loops 2 and 3 are parameterized to satisfy the coordinate angle at 
the boundary points, the grid is constructed as follows. If the grid is to 
contain M pseudo-radial lines (such as line G'G of Fig. 1) and N pseudo-radial 
lines (such as line JRS), the values of the pseudo-radial coordinate are 

r (i ) = i / ( ISI - I ) i = 0, I, 2, 

and the values of the pseudo-azimuthal coordinate are 

s(j) = j / ( M - I) j = 0,1, 2, M-l 

Then the position vector; i.e., the grid locations x, y for each point in the 
grid is given by Eq. (4). 

The preceding has assumed a uniform spacing in the radial direction. 

If radial grid point concentration is desired, it is simply necessary to assume 
a radial distribution function. The present analysis assumed a distribution 
function 

f tonhD(l-r) 1 

R, ['~ lonhD -J «> 

which concentrates points in the wall region. Grid points are then chosen at 
r(i) = (i)/(N-l) and the analysis proceeds as outlined. 
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Grids constructed with this technique have been used successfully in 
conjunction with a Navier-Stokes calculation procedure for a variety of airfoil 
calculations. Typical results are presented in Refs. 5-9. 

User Instruction 

Based upon the input which is discussed in detail in the subsequent 
section, the code constructs a f C'-type coordinate system. In regard to the 
input, the most critical item is the airfoil definition. The airfoil can be 
defined in SUBROUTINE PROFIL which currently can be used to define NACA 4-digit 
airfoils with mean line parameters XCMX and YCMX, and thickness parameter T 
as given by Abbott and Von Doenhoff in Sections 6.3 and 6.4 (Ref. 10). 

The subroutine prints the airfoil type upon the first call to the subroutine at 
FORMAT 1000. Other airfoils in this class can be obtained by varying XCMX, 

YCMX and T. In all cases, the FORMAT statement should be changed by the user 
to reflect the relevant airfoil. 

The subroutine can be modified for other airfoil types. In all cases X is 
input to the subroutine and gives a dimensionless chordwise location, 

0 < x < 1. If NAIR =1, the lower surface coordinate normalized by chord, XAIR 
and YAIR, must be calculated. If NAIR = 0, the upper surface coordinates must 
be calculated. Alternatively, the coordinates can be input as discrete points 
via a NAMELIST AUTO. In these cases care must be taken to adequately define 
the airfoil desired. Presmoothing of the points by user-supplied smoothing 
routine may be beneficial. 

In regard to output, the code prints metric data on a loop-by-loop basis; 
i.e., at all grid points along a loop such as JRS. In assessing the output, 
the critical item is smoothness of variations of the angle between coordinate 
lines and of the Jacobian. Large changes in either may indicate an unsuitable 
grid. This is particularly true of the coordinate angle. In addition, if the 
coordinate angle becomes less than 45° or greater than 135°, the grid may be 
undesirable for computation. Finally, care must be taken to obtain adequate 
resolution. Obviously, this will depend upon the problem considered. However, 
in general the number of pseudo-azimuthal points should not be less than 81 and 
the number of pseudo-radial points not less than 30 for viscous flows. Most 
calculations run to date have utilized grids of 141 x 39 points. 
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Implementation Instruction 

The coordinate generation code has been run successfuly on the NASA 
Langley Research Center CYBER 203 machine using a VS0S2-1 operating system. 

A listing of sample control cards is found at the end of this subsection. 

The first three cards represent the accounting control cards. The fourth card 
attaches the old program library. This is followed by the UPDATE directive 
which creates a compile file. The next card, the TOVPS card, transfers files 
to the 203 for execution. 

This set is followed by a 7/8/9 card indicating end of record, the UPDATE 
directives and another 7/8/9 card. This second 7/8/9 card is followed by the 
CYBER 203 control cards. The first two give accounting information. This is 
followed by the ATTACH which attaches the COMPILE file and the FORTRAN card. 
This is followed by a tape request, a LOAD and a GO. The next two cards copy 
the generated geometry file GE0M2 into GE0012E and then saves GE0012E as a 
permanent file. 


ACCOUNTING CARD 1 
ACCOUNTING CARD 2 
ACCOUNTING CARD 3 
GET(0LDPL = CRDSRA) 

UPDATE, P,C,F,L = A1234. 

TOVPS ( INPUT, C6UD=COMPILE ,UN=****** ,PW=******,AC=******) 
7/8/9 

*C COORD, INTEP, PRFL 
7/8/9 

ACCOUNTING CARD 4 
ACCOUNTING CARD 5 

RESOURCE (TL=50 ,WS=57 6 ,LP=4 , JCAT=SMBAT ) 

ATTACH, COMPILE. 

FORTRAN (L=0UTPUT/ 1000 ,I=C0MPILE ,B=BINARY ,0PT= BL) 
REQUEST(TAPE9/200,RT=W) 

LOAD (B INARY, CN = GO/200 ,L=0UTPUT,GRLPALL= ) 

GO. 

C0PY(GE0M2 ,GE0012E) 

DEFINE, GE0012E. 

7/8/9 

&AUT0 

C4 = 3. , DAMP=6.8 , NRAD = 39, NTH ETA = 141, RMAX0 = 3., 
RMAX1 = .5, JDATA = 1, CCUSP = 0.9, TS = 6., 

IFLAG = 1, NE = 73, LED = 37, 

ICRWR = 1, 8*0, 1, 9*0, 1, 9*0, 1, 4*0, 1, 3*0, 1, 

Rl = 1.0E-06 
SEND 
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Input 


The input is found in NAMELIST $AUTO. It consists of geometry 
specification, grid specification and print flags. The items are 


Geometry Specification 


IFLAG 


DMX(I), DMY(I) 


NE 

LED 

C4 


CCUSP- 


1 Airfoil shape specified via analytic 
formulation in SUBROUTINE PRFL 

2 Airfoil shape specified via discrete data 
points DMX, DMY 

x and y coordinates for discrete data point 
specification. 1=1, NE. The leading edge 
is at I = LED* Begin specification at 
trailing edge x = 1, and proceed in a 
counterclockwise manner over upper surface, 
through I = LED, over lower surface and back 
to trailing edge. Only required when 
IFLAG = 2 

See DMX, DMY 
See DMX, DMY 

Length of branch cut - recommended 
value =3.0 

Location where cusped trailing edge begins - 
recommended value =0.90 
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Grid Specification 


NRAD 

NTHETA 

RMAXO 


RMAX1 

TS 


XLE, MLE 


R1 


Number of pseudo-radial points 

Number of pseudo-azimuthal points 

Distance from branch cut to outer loop 
(h 3 in Fig. 2) 

Secondary loop parameter - set as RMAXO/6 

Parameter to control leading edge 
pseudo-azimuthal resolution. Increasing TS 
increases leading edge resolution relative 
to resolution elsewhere. Recommended 
value = 6.0 

Parameters to set point of maximum leading 
edge resolution 

MLE =2 - Maximum resolution at X = 0 

MLE - 1 - Maximum resolution at X = XLE 
on upper surface 

MLE =0 - Maximum resolution at X = XLE 

on lower surface 

Approximate location of first point from the 
airfoil = Rl*RMAX0 


DAMP - If Rl is not set in NAMELIST $AUT0 

pseudo-radial stretch is set by DAMP. 
Increasing DAMP increases resolution near 
airfoil surface. Recommendation is to set 
Rl 

ANGROT - Allows rotation of inner loop relative to 

outer loop by ANGROT radians; i.e., in 
Fig. 1 loop A f ED f will rotate by 
ANGROT radians 
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Print Flags 


ICRWR(I) 


JDATA 


= 1 Write metric data for loop I on ouput file 

“ 0 Do not write metric data for loop I on 

output file 

= 1 Write inner and outer loop data on output 

file 

= 0 Do not write inner and outer loop data on 

output file 


A sample set of input is given on page 7. 


Output 


The first output item of the code defines the airfoil. The airfoil is 
defined by name if an analytic formulation is used or by discrete coordinates 
if a set of discrete specifying data is used. This is followed by 
NAMELIST $AUTO and, if Rl is set, by the iteration to determine the 
pseudo-radial stretching parameter. This is followed by the location of the 
centerpoint used in constructing the pseudo-circular frontal coordinate 
system. If JDATA = 1, the x and y coordinates as well as the parameter, 
associated with each point defining the inner and outer loops follow. 

If JDATA = 0, this is omitted. 

The next items appearing are the x and y coordinates of each of the four 
loops used to construct the mesh. ARC1 and ARC3 are the inner and outer 
loops, respectively. ARC2 and ARC4 are the secondary or intermediate loops. 
This is followed by the radial distribution function which gives the 
dimensionless radial distance as a function of grid point. Finally, metric 
data are presented for each loop, I, for which ICRWR(I) = 1. The items printed 
are the x and y coordinates (PX and PY), the components of the metric tensor 

Srr» Stt* Srt» the Jacobian of the transformation, and the angle between 
coordinate lines where 
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dx dx dy dy 

d£ drj d£ dr) 


The last four columns are the derivatives 

dx/dr), dy/drj , dx/d£ , dy/d£ 

The final output item is written on TAPE9 and stores the mesh 
coordinates. These are written one loop at a time as xj, y^, i = 1, NTHETA 
as a binary write. 
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Environmental Characteristics 


The program was developed on the CYBER 170, and various versions have been 
used on both the CRAY1 and CYBER 203. A typical grid of 140 x 39 points 
requires less than 10 system time units on the CYBER 203. 
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Figure 1 - Sketch of constructive coordinate system 
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